Elimizde 100m lik tel olduğunu varsayalım. Bu teli bahçenin çevresini kapatmakta kullanmak istiyoruz. Bu şekilde oluşturacağımız bahçelerden hangisinin alanı en büyük olur?
Bu problemde olduğu gibi bir çok farklı türdeki problem optimizasyon yapmamızı gerekli kılar. Örneğin bir fonksiyonun köklerini bulmak isteyebiliriz. Ya da o fonksiyonu minimum ya da maksimum yapan değerleri.
Bir beslenme uzmanı, iki farklı gıda türünü birleştirerek içeriğinde en az 8 birim A vitamini ve 10 birim C vitamini bulunduracak şekilde bir karışım elde etmek istiyor.
Birinci gıda türünde kilogram başına 2 birim A vitamini ve 1 birim C vitamini varken, ikincisinde kilogram başına 1 birim A vitamini ve 2 birim C vitamini vardır.
Birinci gıdanın kilogram maliyeti 50 lira, ikincisinin kilogram maliyeti ise 70 liradır. Karışımın maliyetini mümkün olduğunca düşük tutmak için, bu gıdalardan ne kadar almalıdır?
R’da çok fazla sayıda optimizasyon işi yapan paket vardır. CRAN taskview Optimization sayfasında bunların bir özeti sunulmuştur.
https://cran.r-project.org/web/views/Optimization.html
Bu paketler farklı farklı amaçlar için hazırlanmıştır. Farklı optimizasyon problemleri için farklı paketleri kullanmak gerekir.
Base R daki optimizasyon fonksiyonları saymak istersek
fonksiyonlarını söyleyebiliriz.
Elimizde 100m lik tel olduğunu varsayalım. Bu tel ile bir bahçenin çevresini kapatmak istiyoruz. Bu şekilde oluşturacağımız bahçelerden hangisinin alanı en büyük olur?
Dikdörtgen şeklindeki bahçenin kenarları x,y uzunluğunda olsun. Çevresi \(2(x+y)=100\) olacağından \(y=50-x\) bulunur.
Alanı ise \(xy=x(50-x)=50x-x^2\) şeklindedir. Bu durumda alan fonksiyonu olan \(50x-x^2\) yi maximum yapan x değerini bulmamız gerekmektedir.
x hakkında bildiğimiz şeyler \(x>0\) ve \(x<100\) olduğudur.
Bunlara kısıt (constraint), alan fonksiyonuna ise hedef fonksiyon (objective function) denir.
\(0-100\) arasında alan fonksiyonunun grafiğini çizerek en yüksek olduğu değerin \(x=25\) olduğunu görebiliriz.
curve(50*x-x^2, from = 0, to = 100) abline(v=25)
Alan fonksiyonunu optimize etmek için optimize fonksiyonunu kullanabiliriz.
alan <- function(x) 50*(x)-x^2 opt<-optimize(alan,c(0,100),maximum=TRUE) opt
## $maximum ## [1] 25 ## ## $objective ## [1] 625
Hedef fonksiyonu maksimum yapan değer \(x=25\) ve hedef fonksiyonun bu aralıktaki maksimum değerinin \(625\) olduğu görülür.
Ya da optim fonksiyonunu kullanarak aynı sonuca ulaşabiliriz. Aradaki fark optim fonksiyonunun çok değişkenli optimizasyon da yapabilmesidir.
alan <- function(x) 50*(x)-x^2 opt<-optim(0,alan,method = "CG",control=list(fnscale=-1)) opt
## $par ## [1] 25 ## ## $value ## [1] 625 ## ## $counts ## function gradient ## 90 31 ## ## $convergence ## [1] 0 ## ## $message ## NULL
Dikdörtgenler prizması şeklindeki bahçenin kenarları x,y,z uzunluğunda olsun. Çevresi \(4(x+y+z)=100\) olacağından \(x+y+z=25\) ve \(z=25-x-y\) bulunur.
Bu durumda hacim fonksiyonu olan \(xyz=xy(25-x-y)\) yi maximum yapan x ve y değerlerini bulmamız gerekir.
İki değişken olduğu için optim kullanmalıyız.
hacim <- function(x)
x[1]*x[2]*(25-x[1]-x[2])
opt<-optim(c(0,1),hacim,
control = list(fnscale = -1))
opt$par
## [1] 8.333333 8.333333
opt$value
## [1] 578.7037
Görüldüğü gibi optimum değerlerin x=8.33, y=8.33 ve z=8.33 olduğu bulunur. Yani aslında hacimin maksimize olduğu durum bir küptür.
surf <- function(x, y){x*y*(25-x-y)}
x <- y <- seq(0, 100, length = 50)
z <- outer(x, y, surf)
plotly::plot_ly(x = x, y = x, z = z) |> plotly::add_surface()
optim fonksiyonunun parametreleri şu şekildedir.
optim(par, fn, gr = NULL, ...,
method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent"),
lower = -Inf, upper = Inf,
control = list(), hessian = FALSE)
Aynı problemi 3. bir yoldan çözelim. nlm fonksiyonu Newton metodunu kullanarak minimizasyon yapar.
hacim <- function(x) -1*x[1]*x[2]*(25-x[1]-x[2]) opt<-nlm(hacim,c(0,1)) opt
## $minimum ## [1] -578.7037 ## ## $estimate ## [1] 8.33333 8.33333 ## ## $gradient ## [1] -1.177341e-05 -8.335522e-06 ## ## $code ## [1] 1 ## ## $iterations ## [1] 10
Aynı problemi 4. bir yoldan çözelim. Hacim fonksiyonunu bulurken \(4(x+y+z)=100\) eşitliğini kullanarak z yerine \(25-x-y\) yazmış ve \(xyz=xy(25-x-y)\) bulmuştuk. Bu şekilde 1 değişkeni optimizasyon dışı bırakmıştık.
z’yi devre dışı bırakmadan problemi çözmek için kısıt (constraint) kullanmamız gerekir. Bu durumda hacim fonksiyonu \(xyz\) yi maximum yapmalıyız öyle ki \(4(x+y+z)=100\) olsun. İşte buradaki son eşitlik \(4(x+y+z)=100\) bir kısıt ifadesidir.
Yukarıda bahsedilen fonksiyonlar constrained optimization yapamazlar. Bu şekilde bir optimizasyon için base R’da constrOptim fonksiyonu kullanılabilir.
Bizim kullanmak istediğimiz kısıt \(4(x+y+z)=100\) ya da \(x+y+z=25\) kısıtı bir eşitlik kısıtı(equality constraint) dır. \(x+y+z<100\) gibi bir kısıt olsaydı buna eşitsizlik kısıtı diyecektik.
Sorun şu ki 
constrOptim kullanabilmek için problemimizi değiştirelim. Yine hacmi \(xyz\) maksimize etmek isteyelim, ama bu sefer kısıtımız \(x+y\le{z}\) olsun. Ayrıca problemin feasible olabilmesi için \(0<x,y,z<30\) olma kısıtlarını da ekleyelim. Bunu aşağıdaki gibi yapabiliriz.
hacim <- function(x) x[1]*x[2]*x[3] opt<-constrOptim(f = hacim,theta = c(0,0.1,0.1),grad=NULL, control = list(fnscale = -1), ui=rbind(c(-1,-1,1),c(-1,0,0),c(0,-1,0),c(0,0,-1),c(1,0,0),c(0,1,0),c(0,0,1)), ci = c(0,-30,-30,-30,0,0,0) - 1e-6) opt
## $par ## [1] 14.99994 15.00006 30.00000 ## ## $value ## [1] 6750.001 ## ## $counts ## function gradient ## 950 NA ## ## $convergence ## [1] 0 ## ## $message ## NULL ## ## $outer.iterations ## [1] 3 ## ## $barrier.value ## [1] 0.0264519
Kısıtılarımız \(x+y\le{z}\) ve \(0<x,y,z<30\) şeklinde basit kısıtlar olmasına rağmen neden bu kadar karmaşık ifadeler girdik diye düşünebilirsiniz.
ui=rbind(c(-1,-1,1),c(-1,0,0),c(0,-1,0),c(0,0,-1),c(1,0,0),c(0,1,0),c(0,0,1))
ifadesi aslında \(-x-y+z\ge{0}\) ve \(-x>-30\),\(-y>-30\),\(-z>-30\),\(x>0\),\(y>0\),\(z>0\) eşitsizliklerine karşılık gelir.
ci = c(0,-30,-30,-30,0,0,0) - 1e-6
ifadesi de eşitsizliklerin sağ tarafına aittir. -1e-6 çok küçük bir değeri çıkartarak algoritmanın çalışmasını sağlıyoruz. Başlangıç değerlerini seçmek de önemlidir. \(ui*theta-ci\ge{0}\) olacak şekilde seçilmelidir.
theta = c(0,0.1,0.1)
Base R’daki bu paketler artık biraz eskimiştir. Bunların yerine modern paketler kullanmak daha uygundur. Örnek vermek gerekirse
Dahası için CRAN taskview e bakılabilir.
MASS paketindeki Pima.tr verisinde 21 yaş üzerindeki Pima Indian kadınlarındaki glukoz, bmi, yaş, deri kıvrım kalınlığı, hamilelik sayısı gibi veriler vardır. Bu verisetini kullanarak glukoz değerlerini bmi,yaş ve skin değerlerinin bir fonksiyonu
\(glu=\beta_0+\beta_1*age+\beta_2*bmi+\beta_3*skin\)
olarak yazmamızı sağlayan beta katsayılarını bulmak istiyoruz.
Bu bir regresyondur ve aslında regresyon da bir optimizasyon problemidir. Hedef fonksiyonumuz artıkların kareler toplamıdır. Amacımız bu kareler toplamını minimize etmektir (least squares estimates of regression parameters).
Optimx paketi ile bu işlemi yapalım. Aslında burada optimx fonksiyonuna özel bir işlem kullanmadık. Aynı kodu optim ile de çalıştırabiliriz.
data(Pima.tr,package = "MASS")
min_residuals <- function(data,par) {
with(data, sum((par[1]+par[2]*age+par[3]*bmi+par[4]*skin-glu)^2))}
optimx::optimx(par=c(0,0,0,0), fn=min_residuals, data=Pima.tr)
## p1 p2 p3 p4 value fevals gevals niter ## Nelder-Mead 66.85973 0.9050753 0.7697303 0.1086440 169876.1 459 NA NA ## BFGS 66.90398 0.9047066 0.7678350 0.1097773 169876.1 46 7 NA ## convcode kkt1 kkt2 xtime ## Nelder-Mead 0 TRUE TRUE 0.007 ## BFGS 0 TRUE TRUE 0.001
Optimx sonuçlarında kullanılan algoritmalar (Nelder-Mead ve BFGS) ve elde edilen sonuçlar görülmektedir. p1-p4 değerleri bizim aradığımız katsayılardır.
## p1 p2 p3 p4 value fevals gevals niter ## Nelder-Mead 66.85973 0.9050753 0.7697303 0.1086440 169876.1 459 NA NA ## BFGS 66.90398 0.9047066 0.7678350 0.1097773 169876.1 46 7 NA ## convcode kkt1 kkt2 xtime ## Nelder-Mead 0 TRUE TRUE 0.007 ## BFGS 0 TRUE TRUE 0.000
Bu katsayıları doğrusal regresyon komutu sonuçları ile karşılaştıralım.
Aşağıda doğrusal regresyon komutu sonuçları görülmektedir. Karşılaştırdığımızda katsayıların doğrusal regresyon komutu sonuçları ile çok benzer olduğunu görürüz.
## ## Call: ## lm(formula = glu ~ age + bmi + skin, data = Pima.tr) ## ## Residuals: ## Min 1Q Median 3Q Max ## -77.008 -17.497 -2.838 18.654 82.740 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 66.9040 12.7382 5.252 3.90e-07 *** ## age 0.9047 0.1967 4.599 7.59e-06 *** ## bmi 0.7678 0.4531 1.694 0.0918 . ## skin 0.1098 0.2427 0.452 0.6515 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 29.44 on 196 degrees of freedom ## Multiple R-squared: 0.1487, Adjusted R-squared: 0.1357 ## F-statistic: 11.42 on 3 and 196 DF, p-value: 6.22e-07
Constrained optimizasyonun ne olduğu görmüştük.
constrOptim paketi doğrusal eşitsizlik kısıtlarını kabul etmektedir. Yani \(ax+by\ge{cz}\) şeklindeki kısıtlar.
Eğer kısitlarınız \(ax^2+by^2\ge{cz^2}\) şeklinde ise bu tür doğrusal olmayan kısıtları (nonlinear constraints) kullanabileceğiniz optimizer lar da vardır.
Örnek olarak alabama, nloptr, NlcOptim paketleri söylenebilir.
Himmelbau fonksiyonunun \(f(x,y)=(x^2+y-11)^2+(x+y^2-7)^2\) köklerini bulmak isteyelim.
Başlangıç parametresi \(x_0=(1,1)\) e karşılık gelen kök=(3,2)
## ## Call: ## nloptr::nloptr(x0 = x0, eval_f = himmelbau, eval_grad_f = grad, ## opts = opts) ## ## ## Minimization using NLopt version 2.7.1 ## ## NLopt solver status: 1 ( NLOPT_SUCCESS: Generic success return value. ) ## ## Number of Iterations....: 13 ## Termination conditions: xtol_rel: 1e-08 ## Number of inequality constraints: 0 ## Number of equality constraints: 0 ## Optimal value of objective function: 1.87779471691287e-23 ## Optimal value of controls: 3 2
Başlangıç parametresi \(x_0=(1,-1)\) e karşılık gelen kök=(3.58,-1.84)
## ## Call: ## nloptr::nloptr(x0 = x0, eval_f = himmelbau, eval_grad_f = grad, ## opts = opts) ## ## ## Minimization using NLopt version 2.7.1 ## ## NLopt solver status: 1 ( NLOPT_SUCCESS: Generic success return value. ) ## ## Number of Iterations....: 16 ## Termination conditions: xtol_rel: 1e-08 ## Number of inequality constraints: 0 ## Number of equality constraints: 0 ## Optimal value of objective function: 9.34375355312219e-23 ## Optimal value of controls: 3.584428 -1.848127
Başlangıç parametresi \(x_0=(-1,1)\) e karşılık gelen kök=(-2.80,3.13)
## ## Call: ## nloptr::nloptr(x0 = x0, eval_f = himmelbau, eval_grad_f = grad, ## opts = opts) ## ## ## Minimization using NLopt version 2.7.1 ## ## NLopt solver status: 1 ( NLOPT_SUCCESS: Generic success return value. ) ## ## Number of Iterations....: 15 ## Termination conditions: xtol_rel: 1e-08 ## Number of inequality constraints: 0 ## Number of equality constraints: 0 ## Optimal value of objective function: 1.93938922882802e-21 ## Optimal value of controls: -2.805118 3.131313
Başlangıç parametresi \(x_0=(-1,-1)\) e karşılık gelen kök=(-3.77,-3.28)
## ## Call: ## nloptr::nloptr(x0 = x0, eval_f = himmelbau, eval_grad_f = grad, ## opts = opts) ## ## ## Minimization using NLopt version 2.7.1 ## ## NLopt solver status: 1 ( NLOPT_SUCCESS: Generic success return value. ) ## ## Number of Iterations....: 15 ## Termination conditions: xtol_rel: 1e-08 ## Number of inequality constraints: 0 ## Number of equality constraints: 0 ## Optimal value of objective function: 1.37552028901818e-22 ## Optimal value of controls: -3.77931 -3.283186
Optimx paketi ve içerisindeki optimx fonksiyonu constrained optimizasyon yapamadım. İnternette yaptığını söyleyenler var olsa da benim örneklerimde olmadı. Bu arada constrOptim paketinin dökümantasyonu da çok sınırlı. Nasıl yapıldığını anlamak için internet kaynaklarını iyice araştırmam gerekti.
Base R paketleri optimizasyon için kullanılabilir ama biraz eski oldukları görülmektedir. optimx daha yeni ve kullanımının daha kolay olduğunu söylemektedir. Her paketin farklı bir işlevi ve kabiliyeti olduğu unutulmamalıdır.
